On the origin of non-Gaussian statistics in hydrodynamic turbulence 
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Turbulent flows are notoriously difficult to describe and understand based on first principles. One 
reason is that turbulence contains highly intermittent bursts of vorticity and strain-rate with highly 
non-Gaussian statistics. Quantitatively, intermittency is manifested in highly elongated tails in the 
probability density functions of the velocity increments between pairs of points. A long-standing 
open issue has been to predict the origins of intermittency and non-Gaussian statistics from the 
Navier-Stokes equations. Here we derive, from the Navier-Stokes equations, a simple nonlinear 
dynamical system for the Lagrangian evolution of longitudinal and transverse velocity increments. 
From this system we are able to show that the ubiquitous non-Gaussian tails in turbulence have their 
origin in the inherent self-amplification of longitudinal velocity increments, and cross amplification 
of the transverse velocity increments. 



Intermittency in turbulent flows refers to the vio- 
lent and extreme bursts of vorticity and rates of strain 
that occur interspersed within regions of relatively quiet 
flowQ. These infrequent, but extreme events are be- 
lieved to cause observed deviations from the classical 
Kolmogorov theory of turbulence^- Intermittency also 
has a number of practical consequences since it can lead 
to sudden emergence of strong vortices in geophysical 
flows 0, to modifications of the local propagation speed 
of turbulent flames Q|, etc. One of the observable man- 
ifestations of intermittency is the tendency of velocity 
increments, i.e. the difference between velocities at two 
spatial points separated by a distance £, to display highly 
non-Gaussian statistics when £ is smaller than the flow 
integral scale, L. The tails of velocity-increment prob- 
ability density functions (pdf) are observed to be ex- 
ponential and even stretched exponential^, Q. More- 
over, an inherent asymmetry develops in the distribution 
of the longitudinal velocity increments, i.e. the differ- 
ence of the velocity component in the direction of the 
displacement between the two points. This asymmetry 
yields the well-known negative skewness of longitudinal 
velocity increments 0. While the negativity of skewness 
can be derived from the Navier-Stokes (N-S) equations in 
isotropic turbulence 0, a straightforward mechanistic 
explanation of the origins of stretched exponential tails, 
intermittency, and asymmetry has remained elusive. 

In one dimension for the Burgers equation, the emer- 
gence of negative skewness and long negative tail in the 
pdf starting from random initial conditions is well under- 
stood based on the nonlinear term's tendency to steepen 
the velocity gradient. In 3D turbulent flows, the notion 
of nonlinear "self- amplification" as the cause of intermit- 
tency has long been suspected jfj. Yet, these expectations 
have eluded quantitative analysis due to the difficulty in 
deriving lower-dimensional models that maintain the rel- 
evant information about the vectorial nature of the full 
3D dynamics. Many surrogate models have been pro- 
posed, such as shell models |q, the mapping closure \ 



etc, but the connection with the original N-S equations 
is typically based on qualitative and dimensional resem- 
blances instead of on systematic derivation. 

We consider the coarse-grained N-S equations filtered 
at scale A comparable (and larger) than the scale i. Let 
Ui be the filtered velocity field. Defining the velocity 
gradient tensor Aji — dvii/dxj and taking the gradient 
of the filtered N-S equations one obtains [fl. Hof that the 
rate of change of the velocity gradient is given by 



-4, 



(A jk A kl + 2Q/3 Sji) + Hj 



(1) 



AmnA nm /2 arises from continuity. The 



where Q 
tensor H 
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ji contains the trace-free part of the pressure 
and viscous force gradients 0, Il2j |: 



in which p is the filtered pressure divided by density and 
v the viscosity, = UjUj — uiUj is the subgrid-scale 
(SGS) stress. The time derivative ( ) is a Lagrangian 
material derivative defined as the rate of change of the 
gradient tensor following the local smoothed flow. Set- 
ting Hij = yields the so-called "Restricted Euler" 
dynamics P. H(3| . A fruitful method to model the effects 
of Hij has been to track material deformations usin g ei - 
ther tetrad dynamics^] or the Cauchy-Green tensor|l4| . 
Here we focus on a simpler object - a line element, aiming 
at identifying the mechanism generating intermittency. 
Thus, consider two points separated by a displacement 
vector r of length smaller than, or of the order of, A 
so that the local velocity field is smooth enough to be 
approximated as a linear field. The velocity increment 
between the two points over the displacement r is then 



Sm (r,t) = u t (x + r) - m (x) w Au r k . 



(2) 



The longitudinal and transverse velocity increments, 
8u(r,t) and 6v(r,t) respectively, can be evaluated from 
the two projections of velocity increment Eq. Elinto di- 
rections longitudinal and transverse to r (see FIG.^): 



Su(r,t) = A ki r k 



Sv{r,t) = \P ij {r)A kj r k \ , (3) 
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where Pair) 



TiTj/r 2 and r = |r| 



Note that Su(r,t) and Sv(r, t) correspond to velocity 
increments over a displacement that is evolving, in 
a local linear flow, according to equation fj = A mi r m . 
To study the evolution of velocity increments at a fixed 
scale £, it is necessary to eliminate effects from the chang- 
ing distance between the two points. Consider a line 
that goes through the two points. Still within the as- 
sumption of a locally linear velocity field, the velocity 
increments across a fixed distance £ along this line are 
Su = 5u(r,t)£/r, Sv = Sv(r,t)£/r (see FIG.[TJ). 



Su^u/x+rJ-u/xJ 




Su=Sur.l/i 2 



FIG. 1: Illustrative sketch of velocity increment 5ui(v) be- 
tween two points x(t) and x(t) + r(t), and the components of 
8ui(r) longitudinal and transverse to the displacement vector 
r. The quantities of interest are Su and Sv, defined as the 
components of the velocity increment over a fixed length I. 

Taking time derivatives of Su and Sv, and using the 
expressions for Aji and ri, many terms simplify and one 
arrives at the following "advected delta-vee" system of 
equations: 



Sii = -Su 2 i~ l + Sv 2 l~ x - -Ql + Y, 



Sv 



-2 Su Sv t 



Z. 



(4) 
(5) 



where Y = IH^riTj jr 2 and Z = £HijCjri/r contain the 
anisotropic nonlocal effects of the pressure, inter-scale 
effects of subgrid-scale stresses, and damping effects of 
molecular viscosity (e is a unit vector in the direction of 
the transverse velocity component). The first term on the 
right-hand-side (rhs) of the equation for Sii also occurs in 
ID Burgers equation (the self- amplification effect of neg- 
ative velocity gradients) . The second term indicates that 
the transverse velocity (rotation) tends to counteract the 
self-amplification process. For Sv, the first term on rhs 
of Eq. El suggests exponential growth of Sv at a rate \Su\ 
when Su < 0. This "cross- amplification" mechanism can 
lead to very large values of \Sv\. 

We now pose the question whether the growth of in- 
termittency and the asymmetry of longitudinal velocity 
increments can be understood based on this system of 
equations, but without the effects represented by Y and 
Z (i.e., "Restricted Euler" dynamics). In order to de- 
termine whether this simplified system approximates Sii 



and Sv in real turbulence, comparisons are made with 
direct numerical simulations (DNS) . The rates of change 
of Su and Sv predicted by DNS are obtained by finite 
difference in time from two DNS velocity fields sepa- 
rated by the simulation time-step St = 0.001. The 
data are obtained from a pseudo-spectral simulation of 
the N-S equations, with 256 3 nodes and Taylor-scale 
Reynolds number R\ i=s 162. The velocity fields are 
coarse-grained using a Gaussian filter of characteristic 
length A = 40?7, where r\ is the Kolmogorov length scale, 
yielding filtered velocity fields I*i(x, to) and I*i(x, to + Si) 
(i = 1,2,3). At the initial time to, to every grid-point 
x(to) on the computational mesh, we associate a part- 
ner x(to) + r (^o) at a distance |r(t )| = £ = 40?y in 
some Cartesian direction. For each pair of points we 
measure the longitudinal and transverse velocity incre- 
ments. Then, we find the position to which x(io) and 
x(i ) + r(£o) will be advected by the smoothed velocity 
field, which are, using simple Euler integration, x(*o + 
St) = x(t ) + u(x, to)St, and x(to + St) + r(t + St), where 
r{t + St) = r(t ) + [u(x(to) + r(t ), to) - u(x(t ), to)]** 
is the new displacement vector. The final end-point at a 
fixed distance I is found by moving the material end-point 
x(< + St) + r(to + St) along the new displacement vector 
to the point x.(to + St) + r(to + St)£/\r(to + 5t)\, so that the 
distance is kept fixed. Velocities at the new locations are 
obtained from the stored field at the new time using bi- 
linear interpolation, and the longitudinal and transverse 
components are evaluated, by projections onto direction 
parallel and perpendicular to the new displacement vec- 
tor between the two points. The rate of change of Su 
and Sv is evaluated using first-order finite difference in 
time. Conversely, the rates of change predicted by the 
model system are evaluated as —Su 2 /£ + Sv 2 /£— (2/3)Q£ 
and —2 Su Sv/£ from the measured values of Su, Sv and 
Q. Both real and modeled rates of change are computed 
over a large number of points in the DNS data, and their 
correlation coefficient and joint pdf are evaluated. 
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FIG. 2: Joint pdf of rates of change of velocity increments pre- 
dicted from DNS (filtered at A = AQ-q and taking increments 
over a distance I = 40r;) and the "advected delta-vee" sys- 
tem, (a) longitudinal and (b) transverse velocity increments. 
Results are robust with changes in A and i (with I < A). 

FIG. |3 shows the joint pdf of the model results versus 
the rates of change measured from DNS. A clear correla- 
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tion can be seen between model and DNS results. Corre- 
lation coefficients are 0.54 for the longitudinal and 0.61 
for the transverse velocity increments, indicating that the 
model system captures important (but clearly not all) ef- 
fects seen in the real dynamics. The deviations between 
model system and DNS are caused by the neglected Y 
and Z terms, to be studied in future work. 

After confirming that the simplified system captures 
important trends in 3D fluid turbulence, we explore the 
trends predicted by solutions of the model system. In the 
present Letter we set Q to a constant Qq (numerical tests 
show that allowing Q to evolve in time leads to the same 
short-time behavior to be displayed below, except if Q 
were to be closely correlated with Su and Sv, which is not 
the case in 3D turbulence, since Q depends on velocity 
gradients along two additional directions |17j.). We note 
that for Q — 0, the system describes the relative motion 
of a fluid (with a locally linear velocity field) consisting of 
non-interacting ( "ballistic" ) particles that maintain their 
initial velocity. For Q ^ 0, the particles are subjected 
to a relative force equal to the spherical average of the 
pressure, inter-scale, and viscous damping forces. For the 
case Qo = the analytical solution is 

Su(t) = £e [e t + du £]/{[e t + Su £} 2 + Sv 2 £ 2 }, (6) 
Sv(t) = £ 2 Sv Q e Q /{[e t + Su £} 2 + Sv 2 £ 2 }, (7) 

where eo = Suq + Sv 2 ,. For discrete values of time, this 
defines a mapping (the "advected delta- vee map"). The 
system has an invariant 

U = {5u 2 + 5v 2 )/5v, (8) 

and its (circular) phase-space trajectories are Su 2 + (Sv — 
U /2) 2 = {U /2) 2 , as shown in FIG. El 




FIG. 3: Phase-space portrait of the "advected delta-vee" dy- 
namical system Su — —Su 2 + Sv 2 , Sv = —2 Su Sv (for Qo = 
and £ = 1). 

In order to illustrate the evolution of 6u(t) and Sv(t), 
we start from an ensemble of randomly oriented lines for 
which the velocity increment vectors are initialized from 
a Gaussian distribution. The increments Su(t) and Sv(t) 
over these lines are evaluated at several later times. To 
compare with experimental data, two issues need to be 
considered. First, since Sv(t) is the magnitude of the 



transverse velocity increment vector, it has to be pro- 
jected onto a coordinate direction to obtain a component 
of the transverse increment, Sv c = 6vcos9. For isotropic 
turbulence, the angle 6 between the vector and a fixed 
direction in the transverse plane is uniformly distributed 
in [0, 2n). Therefore, the pdf P£(Sv c ) of Sv c is related to 
that of Sv, P v (8v), by 

1 f + °° dSv 
P v c (Sv c ) = - P v (Sv) (9) 
n J\6v c \ ^JSv^-Sv 2 

Second, an ensemble of randomly oriented lines (with 
uniform measure on a sphere, i.e. a uniform distribu- 
tion of initial solid angles dtto) will tend to concentrate 
along directions of positive elongation. Thus, in order to 
compare model results at later times with data that are 
taken at random directions not correlated with the dy- 
namics, the model results need to be weighted with the 
evolving solid angle measure. Conservation of fluid vol- 
ume implies that £ 3 dflo — r(t) 3 dfl(t), i.e. in directions 
of growing r(t), the solid angle dfl(t) decreases. Thus, 
probabilities must be weighted by 

dn{t)/dn = [£/r{t)} 3 . (io) 

Since f = Sur/£, we can solve for r(t) and then obtain 
dn{t)/dQ. = exp(-3^ _1 f* 5u(t')dt'). Using the solution 
for Su, we obtain df2(*)/dQ = £ 3 [{£ + Su t) 2 + Sv 2 t 2 }- 3 / 2 
for Qo = 0. This factor is used to weight the measured 
time-evolving pdfs from the model system. Note that 
when Svo — ► and Suq < 0, there is an unphysical finite 
time singularity at t — > ^/|<5uo|, when r — ► 0. 

FIG. El shows the evolution of the pdfs of the longitu- 
dinal and transverse velocity increments (both the mag- 
nitude and a component), as time progresses (for the 
case Qq — 0). It is immediately clear that the two main 
qualitative trends observed in turbulence naturally evolve 
from the solution of the system: the skewness towards 
negative values of longitudinal velocity increment, and 
the noticeable flare-up of long tails in the pdfs of trans- 
verse velocity increment. Also, these features appear 
rather quickly: after a non-dimensional time tjr = 0.18 
the pdf is already highly skewed and displays stretched 
exponential tails. Very similar results are observed for 
nonzero values of Qo (using numerical forward time inte- 
gration with a standard fourth order Runge-Kutta rou- 
tine, we tested Qq = ±2): Relative to the results for 
Qo = 0, the pdfs of Su are shifted to the left for Q > 
and to the right for Qo < 0, and only very minor differ- 
ences are seen for Sv. The rapid appearance of stretched 
exponential tails is due to the divergence of the phase- 
space trajectories on the left half of the plane in FIG. El 
For a given initial kinetic energy Suq + Sv 2 , , HSvq is small, 
the invariant Uq can be arbitrarily large. Thus Su and Sv 
can later grow to very large values during the evolution. 

In summary, the model system proves useful in show- 
ing that the emergence of ubiquitous trends of 3D tur- 
bulence, namely intermittent and asymmetric tails in 
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FIG. 4: Evolution of the pdf of velocity increments in time: 

(a) longitudinal velocity increment, (b) the magnitude and 
(c) a component of the transverse velocity increment vector. 
Su is initialized as a standard Gaussian random number and 
Sv as the square root of the sum of the square of two inde- 
pendent Gaussian random numbers, (c) is calculated from 

(b) by numerical integration of Eq. JUJ. For simplicity, we set 
I = 1, so that the characteristic time-scale of the ensemble 
is t — £/5uo\ims = 1. Dotted line in (a) and (c): Gaussian; 
solid: t = 0.03; dashed: t = 0.06; dash-dotted: t = 0.09; 
dash-double-dotted: t = 0.12; long-dashed: t — 0.15; long- 
dash-dotted: t = 0.18; and thin line with squares in (c): 
t — 0.18 without correcting for evolving measure. 



pdfs of velocity increments, occur even in the "ballis- 
tic" case (Qo = 0). Considering all possible random ini- 
tial directions of relative motion, the fraction of particle 
pairs that initially move towards each other is small, thus 
large gradients in small spatial regions occur rather in- 
frequently but are very intense when they occur due to 
the self- amplification mechanism for Su, and the cross- 
amplification mechanism for Sv. While the model system 
thus helps explain the origin and trends towards inter- 



mittency in 3D turbulence, predicting quantitatively the 
level of intermittency remains an open question. It re- 
quires understanding the effects of pressure, inter-scale 
interactions (that depends on interactions of vorticity 
and strains at various scales, see e.g. 0,^3) an d viscos- 
ity that are neglected in the model system. But already, 
the proposed model system could be combined with cas- 
cade, mapping closure, or shell models to enable these 
heuristic approaches to include a more direct link to the 
underlying Navier-Stokes equations. 
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In 2D turbulence, to which the model can also be applied, 
the term —2QI/3 must be replaced with — ^det(A). It 
can be shown that one of the two terms of the determi- 
nant exactly cancels the term — 8u 2 in the equation for 
5u. This cancels the mechanism for growth of negative 
skewness and intermittency in 2D. In 3D there is no full 
cancellation due to the weaker correlations among the 
different directions. More detailed results for the 2D case 
will be reported elsewhere. 



